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Abstract 

In the first part of this paper we give a short review of the hierarchy of 
stochastic models, related to physical chemistry. In the basement of this 
hierarchy there are two models — stochastic chemical kinetics and the 
Kac model for Boltzman equation. Classical chemical kinetics and chem- 
ical thermodynamics are obtained as some scaling limits in the models, 
introduced below. In the second part of this paper we specify some simple 
class of open chemical reaction systems, where one can still prove the exis- 
tence of attracting fixed points. For example, Michaelis - Menten kinetics 
belongs to this class. At the end we present a simplest possible model 
of the biological network. It is a network of networks (of closed chemical 
reaction systems, called compartments), so that the only source of nonre- 
versibility is the matter exchange (transport) with the environment and 
between the compartments. 
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1 Introduction 

Relation between the existing (mathematical) physical theory and future math- 
ematical biology seems to be very intimate. For example, equilibrium is a com- 
mon state in physics but in biology equilibrium means death. Biology should be 
deeply dynamical but this goal seems unreachable in full extent: even in sim- 
plest physical situations the time consuming complexity of any study of local 
dynamics is out of the present state of art. Thus the only possibility would be 
to consider simpler dynamical models (mean field etc.) but to go farther in their 
structure. The obvious first step should have been related to chemical kinetics 
and chemical thermodynamics. Here we present a review of these first results 
and discuss what should be the second step. 

In the first part of this paper we give a short review (in more general terms 
than in [TS]) of the hierarchy of stochastic models, related to physical chemistry. 
In the basement of this hierarchy there are two models — stochastic chemical 
kinetics and the Kac model for Boltzman equation. Classical chemical kinetics 
and chemical thermodynamics are obtained as some scaling limits in the models, 
introduced below. 
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If some physical conditions, as reversibility, are assumed for a closed (without 
matter exchange) system, then we have sufficiently simple behaviour: one can 
prove convergence to a fixed point. However, in many models of physical chem- 
istry and biology, no reversibility condition is assumed, and the behaviour can 
be as complicated as one can imagine. Here we have already some gap between 
physics and biology, and it is necessary to fill in this gap. In the second part 
of this paper we specify some simple class of open chemical reaction systems, 
where one can still prove the existence of attracting fixed points. For example, 
Michaelis - Menten kinetics belongs to this class. At the end we present a sim- 
plest possible model of the biological network. It is a network of networks (of 
closed chemical reaction systems, called compartments), so that the only source 
of nonreversibility is the matter exchange (transport) with the environment and 
between the compartments. 

2 Microdynamics 

Any molecule of mass m can be characterized by translational degrees of freedom 
(velocity v € R 3 , coordinate x G R 3 ) and internal, or chemical (for example, 
rotational and vibrational) degrees of freedom. Internal degrees of freedom 
include the type j — 1, . . . , J of the molecule and internal energy functionals 
Kj(zj),Zj £ Kj, in the space K,- of internal degrees of freedom. It is often 
assumed, see [IT], that the total energy of the molecule i is 



We consider here the simplest choice when Kj is the fixed nonnegative number, 
depending only on j. It can be interpreted as the energy of some chemical 
bonds. 

We consider the set X of countable locally finite configurations X= {xi 7 i>i, ji} 
of particles (molecules) in R 3 , where each particle i has a coordinate Xi , velocity 
vi and type ji . Denote by 93t the system of all probability measures on X with 
the following properties: 

• Coordinates of these particles are distributed as the homogeneous Poisson 
point field of particles on R 3 with some density c. 

• The vectors (vt,ji) are independent of the space coordinates and of the 
other particles. The velocity v of a particle is assumed to be uniformly 
distributed on the sphere with the radius defined by the kinetic energy 
T = m{v\+V2 + v1)/2 of the particle, and the pairs (ji,Ti) are distributed 
via some common density p(J, T) , 



Our first goal will be to define random dynamics on X (or deterministic dy- 
namics on 9Jt). It is defined by a probability space (X '°°,/i), where fi = 



E t =T t + Kj(zj ti ). 
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/i 0,00 is a probability measure on the set X 0,00 of countable arrays X°-°°(t) = 
{xi(t),Vi(t),ji(t)} of trajectories Xi (t) , Vi (t) , ji (t) on intervals Ii = (Ti,T]i), where 
< n < rn < oo. The measure /x belongs to the set of measures 9Jt 0,o ° on 
X°'°°(t), defined by the following properties: 

• If for any fixed < t < oo we denote by p(t) the measure induced by p 
on X, then fj,(t) e OJt. 

• The trajectories Xi(t), Vi (t) ,ji(t) are independent , each of them is a Markov 
process (not necessary time homogeneous). This process is defined by ini- 
tial measure /i(0) on X, by birth and death rates, defining time moments 
Ti,rji, and by transition probabilities at time t, independent of the mo- 
tion of individual particles but depending on the concentration densities 
c t (j, T) at time t. 

• The evolution of the pair (j, T) for the individual particle in-between the 
birth and death moments is defined by the following Kolmogorov equa- 
tions, which control the one-particle process 



defining Markov process with distributions p t (j, T) . The probability kernel 
P depends however on pt(j,T) itself, we shall make it precise below. 

The dynamics we will describe here is based on some earlier mathematical mod- 
els and central dogmas of physical chemistry. The simplest way to rigorously 
introduce the measure p is by the limit of finite volume random dynamics. Ini- 
tial conditions for this dynamics are as follows: at time some number n^ A )(0) 
of molecules are thrown uniformly in the cube A, their parameters (j, T) are 
independent and have some common density p (j,T). not depending on A. Let 
rij(t) — rij A \t) be the number of type j molecules at time t. 

Input-Output (I/O) processes 

Heuristically, our time scale is such that for the unit of time each molecule 
does 0(1) transitions. Then for any macroquantity q of substance its O(q) 
part may change. One should choose time scales for input-output processes 
correspondingly. 

The (output) rate of the jumps rij — > rij — 1 is denoted by , that is with 
this rate a molecule of type j is chosen randomly and deleted from A. Similarly, 
the (input) rate of the jumps rij — > rij + 1 is denoted by Xj\ that is a molecule 
of type j is put uniformly in A with this rate. Dependence of both rates on 
the concentrations can be quite different. To get limiting I/O process after 
(canonical) scaling one can assume that 




A (0) = ^o) Aj A w = ^) A 



(2) 
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where fj,gj > are some functions of all c\ ,Cj , Cj — rij/A are the 
concentrations. However, mostly we restrict ourselves to the case when fj are 
functions of Cj only. In other words, an individual type j molecule leaves the 
volume with rate /j -* (c^ ) / e~- . Denote 

/^if-if. (3) 



Stochastic chemical kinetics 

The hierarchy presented here depends on what parameters of a molecule 
are taken into account. In stochastic chemical kinetics only type is taken into 
account. The state of the system is given by the vector (ni, . . . , nj). There are 
also R reaction types and the reaction of the type r = 1, . . . , R, can be written 

as 

J> jr Af, =0 

3 

where we denote by Mj a j type molecule, and the stoihiometric coefficients 
Vj r > for the products and Vj r < for the substrates. One event of type 
r reaction corresponds to the jump rij — > rij + Vj r , j = 1, . . . , J. Classical 
polynomial expressions (most commonly used) 

A r = A r Y\ n j U3 r 

for the rates of these jumps define a continuous time Markov process, a kind of 
random walk in Zi. This dependence can be heuristically deduced from local 
microdynamics. However, polynomial dependence is not the only possibility, 
see [3T|. Moreover, there can be various scalings for these rates. The scaling 

A r = a r A 7r+1 , 7 r = Vj r 

j:v jr <0 

where a r are some constants and A is some large parameter, is called canonical 
because the classical chemical kinetics equations 

^ = 2>,r(aw) (4) 

r 

for the densities 



Cj (t) = lim A- x 4 A) (t) 

A— foo J 

follow in the large A limit, with some polynomials Rj :T (see below and [1 7J ) . 



Problem 1 It is important to give (at least heuristic) local probabilistic models 
to explain other than polynomial dependence and scalings for the rates. For 
example for arbitrary homogeneous functions as in 
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It is assumed that at time t — as A — » oo, A _1 n( A )(0) — > c(0). 

Chronologically, the first paper in stochastic chemical kinetics was by Leon- 
tovich [21], which appeared from discussions with A.N. Kolmogorov. Other 
references see in [53]. In 70s stochastic chemical kinetics for small R, J was 
studied intensively, see reviews [2j2 , 9J . At the same time the general techniques 
to get limiting equations §Z§ appears in probability theory [321 [5] . Now there are 
many experimental arguments in favor of introducing stochasticity in chemical 
kinetics [TJ[2J[23]. 

Stochastic energy redistribution 

In the classical Kac model [23] the molecules i = l,...,N have the same 
type, but each molecule i has a velocity Vi or kinetic energy T^. In collisions 
the velocities (or the kinetic energies) change somehow. There is still continuing 
activity with deeper results concerning the Kac model, in particular convergence 
rate, see for example [5]. 

One should merge Kac type models with stochastic chemical kinetics. Then 
each molecule i acquires a pair Tj) of parameters: type j and kinetic energy 
T. However this is not sufficient to get energy redistribution. One should intro- 
duce also "chemical" energy. As it is commonly accepted, the general idea is that 
the energy of chemical bonds of a substrate molecule can be redistributed be- 
tween product molecules, part of the energy transforming into heat. To describe 
this phenomena in well-defined terms we introduce fast and slow reactions. Fast 
reactions do not touch chemical energy, that is types, but slow reactions may 
change both kinetic and chemical energies, thus providing energy redistribution 
between heat and chemical energy. 

Examples of reactions: 

1. All chemical reactions are assumed slow — unary (unimolecular) A — >• B, 
binary A + B — > C + D. synthesis A + B -> C, decay C — > A + B etc. In any 
considered reaction the total energy conservation is assumed, that is the sum of 
total energies in the left side is equal to the sum of total energies in the right 
side of the reaction equation. 

2. Fast binary reactions of the type A + B — > A + B, which correspond to 
elastic collisions and draw the system towards equilibrium. 

3. Fast process of heat exchange with the environment, with reactions of the 
type A + B — >• A + B, but where one of the molecules is an outside molecule. 

If there is no input and output, then the Markov jump process is the follow- 
ing. Consider any subset i\ < . . . < i m ( r ) of m{r) = — <o v i* substrate 
molecules for reaction of type r. 

On the time interval (t, t + dt) these molecules have a "collision" with prob- 
ability A'W^-Vbrdt, where b r is some constant. Let the parameters of these 
molecules be jk = j(ik),Tk = T(ik). Denote 

rn rn 
i=l i=l 

and T',K' are defined similarly for the parameters . . . , . . . , T^, of vc\! 



product molecules. The reaction occurs only if 

T + K - K' > (5) 

and then the energy parameters of the product particles at time t + have the 
distribution defined by some conditional density P r (T{, . . . , T' m ,_^\T\, . . . , T m ) 
on the set < T{ + . . . + T' m ,_ x < T + K - K' . By energy conservation then 

m — 1 

T' m , = T + K-K'- T 'i- 

i=l 

This defines a Markov process M^{t) on the finite-dimensional space (note that 
TeR+) 

Q A = (J Rl 1 x ... x fl^ 

(ni,...,n, 7 ) 

where the union is over all vectors (m, . . . ,nj) such that for the array A = 
(Ai, . . . , Aq) of positive integers and for any atom type q = 1, . . . , Q, 

T'I n j a jg = \ 
j 

where a,j g is the number of atoms of type q in the j type molecule. In other 
words, each atom type defines the conservation law A q = const. 

Now, using conditional densities P r , we define the "one-particle" transition 
kernel 

P(t; j\ , Ti \j, T)=Y p(r) ft ii . T i \h T), (6) 

r 

that is the sum of terms P^ corresponding to reactions r, which we define for 
some reaction types. For unimolecular reactions j —> ji the product kinetic 
energy T\ is uniquely defined, thus Pj^j 1 is trivial and for some constants Ujj 1 , 

p(j^ji) = Ujjl 6(T + K — K\ — Ti). 

For binary reactions j, j' — > ji , j[ , 

x c t (i', T') <5(T + if + T' + if' - i - Ti - If ( - T[). 

In particular, for "fast" collisions (which do not change type) we have the same 
transition kernel but with j = = j[. We see that P^' 3 43lJl ' depend on 
the concentrations c t (j, T) = p t (j, T) c(t) . They are defined via the Boltzman 
type equation 

dCtU ^ Tl) =fj{c j ) + Y / {P{t;hm,T)c t {j,T) 

3 

-Pit-J^j^TJctijuT^dT (7) 
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which is similar to the Kolmogorov equation but includes also birth and death 
terms. 

All technicalities about the derivation of the limiting processes see in Ap- 
pendix of [15] . 

Space dynamics 



To get thermodynamics we need also volume, pressure etc. Thus it is neces- 
sary to define space dynamics and also scaling limit. 

In the jump process, defined above, each particle i independently of the 
others, in random time moments 

Ti(u) < hi(u) < . . . < t in (uj) < . . . < <7i(w) 

changes its type and kinetic energy (thus velocity) . For each trajectory ui of the 
jump process we define the local space dynamics as follows. It does not change 
types, energies, velocities, but only coordinates. If at jump moment t of the 
trajectory ui the particle acquires velocity v(uS) = v(t + 0, uS) and has coordinate 
x(t, u>), then at time t + s 

x(t + s, lu) — x(t, lj) + v(lj)s (8) 

unless the next event (jump), concerning this particle, of the trajectory lj occurs 
on the time interval [t, t + s\. We assume periodic boundary conditions or elastic 
reflection from the boundary. We denote this process by X\(t), the state space 
of this process is the sequence of finite arrays X{ = {ji,Xi,Vi}. Thus each 
particle i has a piecewise linear trajectory in the time interval (r,, <7j). 

Theorem 2.1 The thermodynamic limit X '°°(t) — X c (t) of the processes X\(t) 
exists and its distribution belongs to 9J1 0,00 . 

Proof See [T3]. 



3 Scaling limit 

Now we define more restricted (than 9Jt) manifolds of probability measures on 
X: the grand canonical ensemble for a mixture of ideal gases with one important 
difference — fast degrees of freedom are gaussian and slow degrees of freedom 
are constants Kj, depending only on j. 

We consider a finite number nj of particles of types j = 1, . . . , J in a finite 
volume A. Remind that for the ideal gas of the j type particles the grand 
partition function of the Gibbs distribution is 

00 i / nj r r r \ / ™ J m v 2 \ 

9 W)=EdII /. / &i,i<®3,i exp/3 l nj fa- Kj )-Y,-Y*-) 

„ n i' JAJR 3 JI, / \ ■ 1 / 



=0 3 

1 

n 



^ — -(AA) nj exp/3(fj,j — Kj)nj = cxp(AAj exp/3/ij) 



n-j—0 



J ■ 
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where 

A, =/r 3/2 (^) 3/2 , fa = /'., - Kj. 

General mixture distribution of J types is denned by the partition function 
6 = IljLi ©0) P)- The limiting space distribution of type j particles is the 
Poisson distribution with concentration Cj. We will need the formulas relating 
Cj and \if 



Cj a p g ^ Ajexpp^, 



H = /T 1 Inf^AT 1 ) = A*i,o + r 1 Incj + X,, (9) 

where /i^o = — ft luXj is the so called standard chemical potential, it corre- 
sponds to the unit concentration cj = 1. We put c = Ci + . . . + cj. 

We will need Gibbs free energy G and the limiting Gibbs free energy per 
unit volume 

2 = A™ I = ^>^ - 

Define by 9Jt C 9Jt the set of all such measures for any ft, [ii, . . . , and 
by SCto,/3 its subset with fixed /3. 

In the process defined above the kinetic energies are independent but may 
have not x 2 distributions, that is the velocities may not have Maxwell distribu- 
tion. We force them to have it by specifying some trend to equilibrium process 
(elastic collisions) and heat transfer (elastic collisions with outside molecules) 
processes. 

Assume that there is a family M(a), < a < oo, of distributions \i a on R + 
with the following property. Take two i.i.d. random variables £1,^2 with the 
distribution M(a). Then their sum £ = £1 + £ 2 has distribution M(2a). We 
assume also that a is the expectation of the distribution M(a). Denote p(£i|£) 
the conditional density of £1 given £, defined on the interval [0, £]. We put 

p(f)(T 1 \T,T')=p(T 1 \T + T') 

and of course T[ = T + T' — T\. Denote the corresponding generator by . 

We model heat transfer similarly to the fast binary reactions, as random 
"collision" with outside molecules in an infinite bath, which is kept at constant 
inverse temperature ft. The energy of each outside molecule is assumed to have 
X 2 distribution with 3 degrees of freedom and with parameter ft. More exactly, 
for each molecule i there is a Poisson process with some rate h. Denote by 
tik,k = 1,2,..., its jump moments, when it undergoes collisions with outside 
molecules. At this moments the kinetic energy T of the molecule i is trans- 
formed as follows. The new kinetic energy T\ after transformation is chosen 
correspondingly to conditional density p on the interval [0,T + £j/j], where ^ 
are i.i.d. random variables having \ 2 distribution with density cx 1 / 2 exp(— ftx). 
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Denote the corresponding conditional density by P^'(T\\T). In fact, this pro- 
cess amounts to N independent one-particle processes, denote the corresponding 
generator H^K 

Thus we can write the generator as 

H = H(s f ,sp) = H (r) + s f H (f) + s p H {f}] 

where corresponds to slow reactions and Sf,sp are some large scaling fac- 
tors, which eventually will tend to infinity. 

We will force the kinetic energies to become % 2 using the limit Sf — > oo. 

Theorem 3.1 The limits in distribution 

£ c (t) = lim X c (i), O c ,p(t) = lim <t c (t) 

exist for any fixed t. Moreover, the manifold DJIq is invariant with respect to the 
process £ c (i) for any fixed rates u, b, h. The manifolds 9Jto,/3 are invariant with 
respect to D Ci/ g(i). 

Thus, in the process € c (t) the velocities have Maxwell distribution at any 
time moment. For the process O c ^(t) moreover, at any time t the inverse 
temperature is equal to /?, that is there is heat exchange with the environment. 
Our individual molecules still undergo Markov process, but simplified. At the 
same time, the macrovariables undergo deterministic evolution on 9Jlo./3 ■ 

Markov property — chemical kinetics restoration 



Note that initially the jump rates depend on the energies. We show that, 
after the scaling limit, the process restricted on the types will also be Markov. 
We assume that there are only unary and binary reactions but we do not need 
reversibility assumption here. 

Lemma 3.1 The process, projected on types, that is the process (ni(i),..., 
nj{t)) is Markov. It is time homogeneous for unary reaction system and time 
inhomogeneous in general. 

Proof Recall that the jump rates were assumed to have simplest energy 
dependence, that is collisions occur independently of the energies, but reactions 
occur only if energy condition ([5]) is satisfied. Write gp(r) = P(|£| > r) for the 
X 2 random variable £ with inverse temperature j3. 

Assume K\ < . . . < Kj and consider first the case of unary reactions. It 
is easy to see that the process D c> p(t) can be reduced to the Markov chain on 
{1, . . . , J} with rates Vjj> = Ujj> if j > j' , and Vj.y = g^(Ky — Kj)uj.y if j < j' . 
We used here that the kinetic energy distribution is \ 2 a t any time moment. 

Similarly for the binary reaction j,j' — > ji, j[ we define the renormalized 
Markov transition rates as c(j,j' — > ji, j[) — bjjt^.j lj f if Kj + Kj' > Kj 1 +Kj> 
and 

c(j,f h,j[) = b jd ^ jldl P{|6 + 61 > K n + Kj, - (K, + Kj,)} 
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if Kj + Kj> < Kj 1 + Kf . Here are independent and \ 2 with inverse temper- 
ature /?. It is crucial here the use of the scaling limit for fast reactions. 

Thus, in the thermodynamic limit we get the equations without the energies, 
that is the classical chemical kinetics 

^M = ^2R. rm )+f( Cj ). (io) 

r 

Example: monotonicity of Gibbs free energy for closed system with only unary 
reactions 



Assume now that the continuous time Markov chain on {1, . . . , J} with rates 
Ujji is irreducible. We say that this Markov chain is compatible with the equi- 
librium conditions 

m = ... = fjLj (11) 

if its stationary probabilities 7Tj, or stationary concentrations Cj tC = TtjC, satisfy 
the following conditions 

lnci !e + (^1,0 + K x ) = ... = lncxe + (^,/, + Kj). 

Remark This compatibility condition should appear naturally in local dy- 
namics, but it is not clear how to deduce it in the mean field dynamics. Note 
that reversibility is not a sufficient condition for the compatibility condition. 

To exhibit monotonicity for dynamics one needs special Lyapounov functions 
in the space of distributions. For Markov chains this is the Markov entropy with 
respect to stationary measure Ttj, 

7Ti ' 



S M = / J Pj In : 



see for example |14) . 

Recall that the equilibrium function — Gibbs free energy g(t) — undergoes 
deterministic evolution together with the parameters \ij or Cj . We will show that 
at any time moment it coincides with the Markov entropy up to multiplicative 
and additive constants. 

Theorem 3.2 If the compatibility condition ill]) holds, then 

g(t) = »c+-^S M (t) (12) 

and monotone behaviour of the Gibbs free energy density follows. 
Proof We have 

Q ^ ^ ^ 

9 = lim = 2J c iMj = P' 1 Yl C J ln C J + C J (^>° + K i) ( 13 ) 

3 3 3 



/ . ^ " x "0 ^ / j ' 

3 3 



3 



C 
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where the first and the second equalities are the definitions, in the third and 
the fourth equalities we used the formula 

Mj = /T 1 In (^AJ 1 ) = W , + r 1 Inc, + K h (14) 

where 

Ho = -/T 1 In Aj = -/T 1 ( - ^ + In B,-) (15) 

is the so called standard chemical potential, it corresponds to the unit concen- 
tration Cj = 1 for the equilibrium density, see for example [15j. 
At the same time 

S M =J2pjln^ = Cj2cj In — ■ 

We see that for unary reactions one does not need reversibility assumption. 
Monotonicity of Gibbs free energy for closed system with binary reactions 



For binary reactions a similar result holds (we will not formulate it formally) . 
However, we do not have Markov evolution for the concentrations anymore. 
Instead, we have the Boltzman equation for the concentrations, that is the 
so called nonlinear Markov chain on {1, . . . , J}. Then, instead of the Markov 
entropy one should take the Boltzman entropy with respect to some one-point 
distribution p^ (see definitions in |17| ) 




which coincides with the Markov entropy for ordinary Markov chains. For the 
monotonic behaviour of the Boltzman entropy, one should assume reversibility 
or a more general condition — unitarity, called local equilibrium in |17j . Under 
this condition the monotonicity of the Boltzman entropy was proved in |17j . We 
get the same formula as (I12[) if we replace Sm by —Sh- 

Note that under these conditions Pj(t) is a time inhomogeneous Markov 
chain. In fact, in the long run, that is as t — > oo, the transition rates for 
one-particle inhomogeneous Markov chain, in the vicinity of the fixed point, 
is asymptotically homogeneous. This shows that binary case is asymptotically 
close to the unary case. 



4 Open thermodynamic compartments 

Reversible and nonreversible processes 

Our systems in finite volume evolve via Markov dynamics. It is not known 
when and how this dynamics could rigorously be deduced from the local phys- 
ical laws. However, there are many arguments that reversibility is a necessary 
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condition for this. Reversibility is a particular case of the unitarity property 
of the scattering matrix of a collision process. It was called local equilibrium 
condition in [T71 [5] ) • 

The reversibility gives strong corollaries for the scaling limits — 1) Boltz- 
man monotonicity and 2) attractive fixed points. We call chemical networks 
with properties 1) and 2) thermodynamic compartments. Denote the class of 
such systems T. These systems are a little bit more general than the systems, 
corresponding to the systems with local physical laws (in particular, having 
convergence to equilibrium property). For example, any unimolecular reaction 
system belongs to T, because, as we saw above, the Markov entropy is the 
Boltzman entropy here. However, biological systems obviously are not of class 
T. There are different ways to generalize class T systems. 

The first one is quite common: in chemical and biological systems stochas- 
tic processes usually are not assumed to be reversible. However, without the 
reversibility assumption the time evolution could be as complicated as possible 
(periodic orbits, strange attractors etc.). That has advantages — one can adjust 
to real biological situations, and disadvantages — too many parameters, even 
arbitrary functions. Normally, the rate functions Rj^ r can be rather arbitrarily 
chosen, typical example where this methodology is distinctly pronounced is |33| . 
connections with physics lost etc. In other words, theory becomes meaningless 
when one can adjust it to any situation. 

Another way could be a hierarchy of procedures to introduce nonreversibility 
in a more cautious way. Each further step to introduce nonreversibility is as 
simple as possible and each is related to time scaling, for example, reversible 
dynamics is time scaled and projected on a subsystem. We start to study here 
the simplest type of such procedures. In our case the Markov generator will be 
the sum of two terms, 

H H rev -\- H nonTev: (16) 

where the first one is reversible and the other one is not, but the latter cor- 
responds only to input and output processes. One of technical reasons to 
choose such nonreversible hamiltonian is to keep invariance of the manifolds 
2rt,2Ko,9Jto,/3- 

In principle, another philosophy is possible — large deviation or other rare 
event conditioning, this we do not discuss here. 

Example 1: steady states for open unimolecular systems 



We consider the case with J = 2 and unary reactions only, however the 
following assertions help to understand how more general open systems can 
behave. Consider first the thermodynamic limit, and then the stochastic finite 
volume problem. 

In the thermodynamic limit the following equations for the concentrations 
Cj (t),j = 1,2, hold: 

dc! dc 2 

-jr = ~v\c\ + ia>c 2 + h, -jr = ^l c i - ^c 2 + /2, 
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where v\ = U12, 1^2 = "21 an d fj are defined by (J3]). Possible positive (i.e., 
C\,C2 > 0) fixed points satisfy the following system: 

/i(ci) + h{c2) = 0, -v\C\ + v 2 c2 + /i(ci) = 0. 

For example, for constant fj a positive fixed point exists for any c sufficiently 
large and equals 

v 2 c - f 2 vie - fx 
Cl = : , c 2 = ; . 

V\ + V2 V\+ V2 

In the linear case, that is for fj — ajCj, for the existence of a positive fixed point 
it is necessary and sufficient that cij have different signs and \a,j\ < V\-\-v<2- Then 
the positive fixed point is unique and is defined by 

v 2 c 

Cl = ■ . 

v\ + V2 - ai 

For faster than linear growth of fj fixed points cannot exist for large c. 
We see from these formulas that the equilibrium fixed point 

ci = ■ , c 2 = ■ 

V\ +V2 V\ + V 2 

(for the corresponding closed system) is slightly perturbed if fj (or a,j) are 
small. Moreover, the perturbed fixed point is still attractive. This is true in 
more general situations as well. 

Now consider the stochastic (finite volume) case. 

Proposition 4.1 Assume that fj are constants. In a finite volume the process 
is ergodic if^2fj < 0j transient if ^2 fj > and null recurrent if^2fj =0. 

Proof Note that the number of particles is conserved and the number 
of states is finite if there is no I/O, otherwise the Markov chain is count- 
able: a random walk on Z\ = { (711,722) : n\U2 > 0}. There are jumps 
(711,712) — > (ni — l,ri2 + 1) or (711,712) — > (771 + l,ri2 — 1) due to reactions, 
denote their rates i^itii, V2n2 correspondingly. There are also jumps (771,712) - ► 
(771 ± 1, 772), (ni, 712) — > (771,772 ± 1) due to input-output with the parameters 
Oj A and bjA correspondingly. 

Transience and ergodicity can be obtained using Lyapounov function nx + n^ 
and the results from [7] . To prove null recurrence note that for sufficiently large 
c the system should be in the neighbourhood of the fixed point, which exists for 
c sufficiently large. Thus one can also use the same Lyapounov function. 

General conclusion is that only null recurrent case is interesting. However, 
models with constant rates are too naive. It is reasonable that there are reg- 
ulation mechanisms which give more complex dependence of fj on the rates. 
Unfortunately, there is no firm theoretical basis to get exact dependence of 
reaction and I/O rates on the densities. 

Example 2: stochastic Michaelis - Menten kinetics 
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The generator for Michaelis - Menten kinetics is of type p^|) only in some 
approximation. This model has 4 types of molecules: E (enzyme), S (substrate), 
P (product) and ES (substrate-enzyme complex). There are 3 reactions 

E + S -> ES, ES -> E + S, ES E + P 

with the rates &i A _1 n£ns, fc_in,£s, A^n^s correspondingly. We can also fix 
somehow the output rate for P and input rate for S. 

If k% = then, as a zero'th approximation, we have a reversible Markov 
chain. In fact, there are conservation laws 

n E + n ES = m(E), n s + n ES = m(S) 

for some constants m(E), m(S). Thus we will have random walk for one variable, 
say n E Si 011 the interval [0, min (m(-E'), m(S))], with jumps ties — > ties i 1- 
Such random walks are always reversible. The stationary probabilities for this 
random walk are concentrated around the fixed point of the limiting equations 
of the classical kinetics 

= kic s c E - {k-i + k 2 )c E s (17) 

defined by 

ces ■ 



cs 



+ bc s 

for some constants a, b, defined by m(E),m(S). If k 2 > but small compared 
to k\, k-i, then up to the first order in ki we have the P production speed 

dc P c s 
—— = k 2 c E s = k 2 - 



dt a + bcs 

We could also look on this kinetics as on the simple random walk. We have to 
introduce (arbitrarily) output rate for the product P and adjust the input rate 
of S so that the system becomes null-recurrent. In fact, due to the conservation 
law tie + ties — m(E) we have random walk on the half strip {(ns,UEs)} — 
Z+ x (0, m(E)). The null-recurrence condition can be obtained using methods 
of [7], we will not discuss this here. 



5 Network of thermodynamic compartments 

We call thermodynamic compartments, introduced above, networks of rank 1. 
We saw that they have fixed points, and thermodynamics plays the central 
role there. It can be some tightly dependent and/or space localized system of 
chemical reactions. 

Network of rank 2 consists of vertices a — networks of rank 1, and directed 
edges, that is compartments are organized in a directed graph. Directed edge 
from compartment a to compartment a' means that there is a matter flow from 
a to a'. Matter exchange between two compartments suggests some transport 
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mechanism. It is natural that there is a time delay between the moments of 
departure from a and arrival to a' . The simplest probabilistic model could be 
the following. Each j type molecule leaves a for the destination a 1 with rates 
fj,a,a'i similar to defined in @, and after some random time r(j,a,a') arrives 
to a'. Times t(J, a, a') are independent and their distribution depends only on 
j,a,a'. One can imagine that there is an effective distance L(a,a') between 
a and a' and some transportation mechanism, which defines effective speed to 
go through this distance. For example, it can be transport through membrane, 
which can be represented as a layer [0, L] x R 2 of thickness L. During time 
t(j, a, a') the particle is absent from the network, it has left a but has not yet 
arrived to a' . 

Denote by c a j the concentration of type j molecules in the compartment a. 
Limiting equations are 

= " /£(<V(*)) + £ WMt ~ r(j, a, a')) 

a 

^ ^ fj,a',a {pa' (0 ) ^ ^ ^a r ,jrRa r ,r {poi 1 (^) ) 
a r 

where c a — (c a ,i, ■ ■ ■ c a . j), f^ j is the input rate to a from external environment, 

f^j is the output rate from a to the external environment. Note that these 
equations are random due to random delay times r. In the first approximation 
one can consider r constant, however random time delays seem very essential to 
restore randomness on the time scale, higher than microscopic, in the otherwise 
deterministic classical chemical kinetics. 

Note that the above written equations follow from a similar microscopic 
model — we will not formally formulate it, because it is obvious from our 
previous constructions: the corresponding manifold is x ae A$Jl a , where A is the 
set of compartments, 9Jl a is the manifold for the compartment a. 

The following problems and phase transitions can be discussed in the defined 
model on the rigorous basis (in progress): 

1. The method of thermodynamic bounds in the thermodynamic networks, 
defined in |30J. 

2. (Phase transitions due to transport rates.) Normal functioning of the 
network can be close to the system {c a , J;e } of equilibrium fixed points in each 
compartment a. Such situation can be called homeostasis. Homeostatic regula- 
tion — keeping the system close to some system {c a j >e }. If there is no transport, 
then the compartments are independent and the fixed points inside them are 
pure thermodynamic. Under some transport rates the fixed points change in a 
stable way, they smoothly depend on the transport parameters. However, under 
some change of the transport rates, the fixed points may change drastically: the 
system goes to other basin of attraction. 

3. (Phase transition due to time desynchronization.) It is known now that 
even a decease can be a consequence of timing errors. For a network of rank 2, 
having for instance a cyclic topology (this is called circuit in [IS]), assume that 
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the input rates change periodically or randomly in time. The question is: to 
what process the concentrations converge and with what speed ? This time 
behaviour could be the next step in the analysis of the structure of logical 
networks in the sense of |29| . 
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